function [AB, AL, AC, AR, AT, Si, Sh, dx, dy, dv, dc, siga, sct, chi] = ...
    mesh_center_2d_coef( dat, numg, numm, xcm, xfm, ...
    ycm, yfm, mt, src, it, BC, IBSL, IBSR, IBSB, IBST)

% function A = twoDcoef( dat, numg, numm, xcm, xfm, ycm, yfm, mt, src, it, BC,
%    IBSL, IBSR, IBSB, IBST)
%   This function produces the coefficient matrix and source for solution
%   of the two-dimensional diffusion equation.
%
%   The script is long, but the structure repeats and should be easy to
%   follow.
%   
%   Si = inhomogeneous (i.e. fixed) source
%   Sh = nuSig*dV (i.e. homogeneous source)
%   
%   modified on 05/16/2010


% expand the bc's
BCL = BC(1);
BCR = BC(2);
BCB = BC(3);
BCT = BC(4);

% define data components
dc = zeros(numm,numg); 
sr=dc; xi=dc; ab = dc; ns = dc;
sc = zeros(numm,numg,numg);
for m = 1:numm
    for g = 1:numg %g+(m-1)*g
        dc(m,g) = dat(numg*(m-1)+g,1);   % diffusion coefficient
        sr(m,g) = dat(numg*(m-1)+g,2);   % removal cross-section
        ns(m,g) = dat(numg*(m-1)+g,3);   % fission cross-section
        xi(m,g) = dat(numg*(m-1)+g,4);   % fission spectrum
        for gg = 1:numg
            sc(m,g,:) = dat(numg*(m-1)+g,5:end);
        end
        if numg > 1
            ab(m,g) = sr(m,g) - sum(dat(numg*(m-1)+1:numg*m, 4+g ));
        else
            ab(m,g) = sr(m,g);
        end
    end
end

% deal with "abdata"         ab(m,g) = abdat(numg*(m-1)+g,1); % absorption cross-section

% number of fine meshes in x and y coordinates
N = sum(xfm);
M = sum(yfm);
% number of coarse meshes in x and y coordinates
CN = length(xfm);
CM = length(yfm);
% compute dx and dy vectors
dx = zeros(N,1);
dy = zeros(M,1);
% coarse mesh index for fine mesh
cix=zeros(N,1);
ciy=zeros(M,1);
j = 0;
for i = 1:CN
    dx( (j+1):(j+xfm(i)) ) = (xcm(i+1)-xcm(i))/xfm(i);
    cix( (j+1):(j+xfm(i)) )=i;
    j = sum(xfm(1:i));
end
j = 0;
for i = 1:CM
    dy( (j+1):(j+yfm(i)) ) = (ycm(i+1)-ycm(i))/yfm(i);
    ciy( (j+1):(j+yfm(i)) )=i;
    j = sum(yfm(1:i));
end
% volume
dv   = zeros( N*M, 1 );
lenk = N*M;
AC   = zeros( lenk, numg );
AL   = zeros( lenk-1, numg );
AR   = AL;
AB   = zeros( lenk-N, numg );
AT   = AB;
Si   = zeros( N*M, numg );       % fine mesh inhomogeneous source
Sh   = Si;                       % fine mesh homogeneous source
siga = Si;                       % fine mesh absorption cross section
chi  = Si;                       % fine mesh chi-spectrum
sct  = zeros( N*M, numg, numg ); % fine mesh scattering matrix     

for i = 1:N
    for j = 1:M
        k = i+(j-1)*( N );
        dv(k)=dx(i)*dy(j);
    end
end

% Note, this is the right general idea.  For the response matrix approach,
% it will be necessary to use the Legendre polynomials again.
% define boundary sources if required
if BCL == 2
    alphaL = 2*IBSL;%rand(M,numg);%';%ones(M,numg);
end
if BCR == 2
    alphaR = 2*IBSR;%ones(M,numg);
end
if BCB == 2
    alphaB = zeros(N,numg); alphaB(:,1)=1;
end
if BCT == 2
    alphaT = zeros(N,numg); alphaT(:,1)=1;
end

    
for g = 1:numg % group loop
    
    % interior coefficients
    for i = 2:N-1
        for j = 2:M-1

            % materials
            mL = mt( cix(i-1), ciy(j  ) );  
            mR = mt( cix(i+1), ciy(j  ) );  
            mB = mt( cix(i  ), ciy(j-1) );  
            mT = mt( cix(i  ), ciy(j+1) ); 
            mC = mt( cix(i  ), ciy(j  ) );
            
            k = i+(j-1)*( N ); % K matrix row index

            AL( k-1, g )   = -2*dc(mL,g)*dc(mC,g)*dy(j)/(dx(i-1)*dc(mC,g)+dx(i)*dc(mL,g));
            AR( k, g )     = -2*dc(mR,g)*dc(mC,g)*dy(j)/(dx(i+1)*dc(mC,g)+dx(i)*dc(mR,g));
            AB( k-N, g )   = -2*dc(mB,g)*dc(mC,g)*dx(i)/(dy(j-1)*dc(mC,g)+dy(j)*dc(mB,g));
            AT( k, g )     = -2*dc(mT,g)*dc(mC,g)*dx(i)/(dy(j+1)*dc(mC,g)+dy(j)*dc(mT,g));
            AC( k, g )     = dx(i)*dy(j)*sr(mC,g)-(AL(k-1,g)+AR(k,g)+AB(k-N,g)+AT(k,g));
            
            % volume times scattering and fission xsecs for multiplication
            % with phi on the rhs of power iterations
            sct(k,g,:) = dx(i)*dy(j)*sc(mC,g,:);
            if it == 1
                Sh( k, g )  = dx(i)*dy(j)*ns(mC,g);  % put fixed source back  
                chi(k,g)    = xi(mC,g);
            else 
                Si( k, g )  = dx(i)*dy(j)*src(cix(i),ciy(j),g); % fixed, group, volume source
                Sh( k, g )  = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
                chi(k,g)    = xi(mC,g);
            end
            siga(k,g)  = dx(i)*dy(j)*ab(mC,g);              

        end
    end
    
    % left and right edges
    for j = 2:M-1

        i = 1; % LEFT EDGE
        k = i+(j-1)*( N );

        if BCL==1
            beta=1;
        else
            beta=0;
        end
         
        % materials
        mR = mt( cix(i+1), ciy(j  ) );
        mB = mt( cix(i  ), ciy(j-1) );
        mT = mt( cix(i  ), ciy(j+1) );
        mC = mt( cix(i  ), ciy(j  ) );

        AR( k, g )     = -2*dc(mR,g)*dc(mC,g)*dy(j)/(dx(i+1)*dc(mC,g)+dx(i)*dc(mR,g));
        AB( k-N, g )   = -2*dc(mB,g)*dc(mC,g)*dx(i)/(dy(j-1)*dc(mC,g)+dy(j)*dc(mB,g));
        AT( k, g )     = -2*dc(mT,g)*dc(mC,g)*dx(i)/(dy(j+1)*dc(mC,g)+dy(j)*dc(mT,g));
        AC( k, g )     = dx(i)*dy(j)*sr(mC,g)...
                         - (AR(k,g)+AB(k-N,g)+AT(k,g))...
                         + dc(mC,g)*dy(j)*2*(1-beta)/(4*(1+beta)*dc(mC,g)+dx(i)*(1-beta));
        
        if it == 1
        	Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
            chi(k,g)   = xi(mC,g);
        else
            Si( k, g ) = dx(i)*dy(j)*src(cix(i),ciy(j),g); % fixed, group, volume source
            Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
            chi(k,g)   = xi(mC,g);
        end
        sct(k,g,:) = dx(i)*dy(j)*sc(mC,g,:); 
        siga(k,g)  = dx(i)*dy(j)*ab(mC,g);               
        if BCL == 2 % add the incident source
            Si(k,g) = Si(k,g) + 4*dc(mC,g)*dy(j)*alphaL(j,g)/(4*dc(mC,g)+dx(i));
        end
        if g==1                
            dv( k, 1 ) = dx(i)*dy(j);
        end            
        
        i = N; % RIGHT EDGE
        k = i+(j-1)*( N );    
        
        if BCR==1
            beta=1;
        else
            beta=0;
        end
        
        % materials
        mL = mt( cix(i-1), ciy(j  ) );
        mB = mt( cix(i  ), ciy(j-1) );
        mT = mt( cix(i  ), ciy(j+1) );
        mC = mt( cix(i  ), ciy(j  ) );

        AL( k-1, g )   = -2*dc(mL,g)*dc(mC,g)*dy(j)/(dx(i-1)*dc(mC,g)+dx(i)*dc(mL,g));
        AB( k-N, g )   = -2*dc(mB,g)*dc(mC,g)*dx(i)/(dy(j-1)*dc(mC,g)+dy(j)*dc(mB,g));
        AT( k, g )     = -2*dc(mT,g)*dc(mC,g)*dx(i)/(dy(j+1)*dc(mC,g)+dy(j)*dc(mT,g));
        AC( k, g )     = dx(i)*dy(j)*sr(mC,g)...
                         - (AL(k-1,g)+AB(k-N,g)+AT(k,g))...
                         + dc(mC,g)*dy(j)*2*(1-beta)/(4*(1+beta)*dc(mC,g)+dx(i)*(1-beta));
        
        if it == 1
        	Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
            chi(k,g)   = xi(mC,g);
        else
            Si( k, g ) = dx(i)*dy(j)*src(cix(i),ciy(j),g); % fixed, group, volume source
            Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
            chi(k,g)   = xi(mC,g);
        end
        sct(k,g,:) = dx(i)*dy(j)*sc(mC,g,:);  
        siga(k,g)  = dx(i)*dy(j)*ab(mC,g);                      
        if BCR == 2 % add the incident source
            Si(k,g) = Si(k,g) + 4*dc(mC,g)*dy(j)*alphaR(j,g)/(4*dc(mC,g)+dx(i));
        end

    end

    % top and bottom edges 
    for i = 2:N-1
        
        j = 1; % BOTTOM EDGE
        k = i+(j-1)*( N );
        
        if BCB==1
            beta=1;
        else
            beta=0;
        end
        
        % materials
        mL = mt( cix(i-1), ciy(j  ) );
        mR = mt( cix(i+1), ciy(j  ) );
        mT = mt( cix(i  ), ciy(j+1) );
        mC = mt( cix(i  ), ciy(j  ) );
        
        AL( k-1, g )   = -2*dc(mL,g)*dc(mC,g)*dy(j)/(dx(i-1)*dc(mC,g)+dx(i)*dc(mL,g));
        AR( k, g )     = -2*dc(mR,g)*dc(mC,g)*dy(j)/(dx(i+1)*dc(mC,g)+dx(i)*dc(mR,g));
        AT( k, g )     = -2*dc(mT,g)*dc(mC,g)*dx(i)/(dy(j+1)*dc(mC,g)+dy(j)*dc(mT,g));
        AC( k, g )     = dx(i)*dy(j)*sr(mC,g) ...
                         - (AL(k-1,g)+AR(k,g)+AT(k,g)) ...
                         + dc(mC,g)*dx(i)*2*(1-beta)/(4*(1+beta)*dc(mC,g)+dy(j)*(1-beta));
                     
        if it == 1
        	Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
            chi(k,g)   = xi(mC,g);
        else
            Si( k, g ) = dx(i)*dy(j)*src(cix(i),ciy(j),g); % fixed, group, volume source
            Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
            chi(k,g)   = xi(mC,g);
        end
        sct(k,g,:) = dx(i)*dy(j)*sc(mC,g,:);       
        siga(k,g)  = dx(i)*dy(j)*ab(mC,g); 
        if BCB == 2 % add the incident source
            Si(k,g) = Si(k,g) + 4*dc(mC,g)*dx(i)*alphaB(i,g)/(4*dc(mC,g)+dy(j));
        end

 
        j = M; % TOP EDGE
        k = i+(j-1)*( N );
        
        if BCT==1
            beta=1;
        else
            beta=0;
        end
        
        % materials
        mL = mt( cix(i-1), ciy(j  ) );
        mR = mt( cix(i+1), ciy(j  ) );
        mB = mt( cix(i  ), ciy(j-1) );
        mC = mt( cix(i  ), ciy(j  ) );
        
        AL( k-1, g )   = -2*dc(mL,g)*dc(mC,g)*dy(j)/(dx(i-1)*dc(mC,g)+dx(i)*dc(mL,g));
        AR( k, g )     = -2*dc(mR,g)*dc(mC,g)*dy(j)/(dx(i+1)*dc(mC,g)+dx(i)*dc(mR,g));
        AB( k-N, g )   = -2*dc(mB,g)*dc(mC,g)*dx(i)/(dy(j-1)*dc(mC,g)+dy(j)*dc(mB,g));
        AC( k, g )     = dx(i)*dy(j)*sr(mC,g) ...
                         - (AL(k-1,g)+AR(k,g)+AB(k-N,g)) ...
                         + dc(mC,g)*dx(i)*2*(1-beta)/(4*(1+beta)*dc(mC,g)+dy(j)*(1-beta));

        if it == 1
        	Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
            chi(k,g)   = xi(mC,g);
        else
            Si( k, g ) = dx(i)*dy(j)*src(cix(i),ciy(j),g); % fixed, group, volume source
            Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
            chi(k,g)   = xi(mC,g);
        end  
        sct(k,g,:) = dx(i)*dy(j)*sc(mC,g,:);
        siga(k,g)  = dx(i)*dy(j)*ab(mC,g);
        if BCT == 2 % add the incident source
            Si(k,g) = Si(k,g) + 4*dc(mC,g)*dx(i)*alphaT(i,g)/(4*dc(mC,g)+dy(j));
        end


    end
    
    % finally, the corners

    % BOTTOM LEFT ---------------------------------------------------------
    i=1; j=1;
    k = i+(j-1)*( N );
    if BCL==1
        beta1=1;
    else
        beta1=0;
    end
    if BCB==1
        beta2=1;
    else
        beta2=0;
    end

    % materials
    mR = mt( cix(i+1), ciy(j  ) );
    mT = mt( cix(i  ), ciy(j+1) );
    mC = mt( cix(i  ), ciy(j  ) );

    AR( k, g )     = -2*dc(mR,g)*dc(mC,g)*dy(j)/(dx(i+1)*dc(mC,g)+dx(i)*dc(mR,g));
    AT( k, g )     = -2*dc(mT,g)*dc(mC,g)*dx(i)/(dy(j+1)*dc(mC,g)+dy(j)*dc(mT,g));
    AC( k, g )     = dx(i)*dy(j)*sr(mC,g) ...
                     - (AR(k,g)+AT(k,g)) ...
                     + dc(mC,g)*dy(j)*2*(1-beta1)/(4*(1+beta1)*dc(mC,g)+dx(i)*(1-beta1)) ...
                     + dc(mC,g)*dx(i)*2*(1-beta2)/(4*(1+beta2)*dc(mC,g)+dy(j)*(1-beta2));

    if it == 1
    	Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
    	chi(k,g)   = xi(mC,g);
    else
    	Si( k, g ) = dx(i)*dy(j)*src(cix(i),ciy(j),g); % fixed, group, volume source
    	Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
    	chi(k,g)   = xi(mC,g);
    end
    sct(k,g,:) = dx(i)*dy(j)*sc(mC,g,:);
    siga(k,g)  = dx(i)*dy(j)*ab(mC,g);
    if BCL == 2 % add the (left) incident source
        Si(k,g) = Si(k,g) + 4*dc(mC,g)*dy(j)*alphaL(j,g)/(4*dc(mC,g)+dx(i));
    end
    if BCB == 2 % add the (bottom) incident source
        Si(k,g) = Si(k,g) + 4*dc(mC,g)*dx(i)*alphaB(i,g)/(4*dc(mC,g)+dy(j));
    end    

    % TOP LEFT ------------------------------------------------------------
    i=1; j=M; 
    k = i+(j-1)*( N );
    if BCL==1
        beta1=1;
    else
        beta1=0;
    end
    if BCT==1
        beta2=1;
    else
        beta2=0;
    end
    
    % materials
    mR = mt( cix(i+1), ciy(j  ) );
    mB = mt( cix(i  ), ciy(j-1) );
    mC = mt( cix(i  ), ciy(j  ) );

    AR( k, g )     = -2*dc(mR,g)*dc(mC,g)*dy(j)/(dx(i+1)*dc(mC,g)+dx(i)*dc(mR,g));
    AB( k-N, g )   = -2*dc(mB,g)*dc(mC,g)*dx(i)/(dy(j-1)*dc(mC,g)+dy(j)*dc(mB,g));
    AC( k, g )     = dx(i)*dy(j)*sr(mC,g) ...
                     - (AR(k,g)+AB(k-N,g)) ...
                     + dc(mC,g)*dy(j)*2*(1-beta1)/(4*(1+beta1)*dc(mC,g)+dx(i)*(1-beta1)) ...
                     + dc(mC,g)*dx(i)*2*(1-beta2)/(4*(1+beta2)*dc(mC,g)+dy(j)*(1-beta2));
                 
    if it == 1
    	Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
    	chi(k,g)   = xi(mC,g);
    else
    	Si( k, g ) = dx(i)*dy(j)*src(cix(i),ciy(j),g); % fixed, group, volume source
    	Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
    	chi(k,g)   = xi(mC,g);
    end
    sct(k,g,:) = dx(i)*dy(j)*sc(mC,g,:);
    siga(k,g)  = dx(i)*dy(j)*ab(mC,g);
    if BCL == 2 % add the left incident source
        Si(k,g) = Si(k,g) + 4*dc(mC,g)*dy(j)*alphaL(j,g)/(4*dc(mC,g)+dx(i));
    end
    if BCT == 2 % add the top incident source
        Si(k,g) = Si(k,g) + 4*dc(mC,g)*dx(i)*alphaT(i,g)/(4*dc(mC,g)+dy(j));
    end    

    % BOTTOM RIGHT --------------------------------------------------------
    i=N; j=1; 
    k = i+(j-1)*( N );
    if BCR==1
        beta1=1;
    else
        beta1=0;
    end
    if BCB==1
        beta2=1;
    else
        beta2=0;
    end
    
    % materials
    mL = mt( cix(i-1), ciy(j  ) );
    mT = mt( cix(i  ), ciy(j+1) );
    mC = mt( cix(i  ), ciy(j  ) );

    AL( k-1, g )   = -2*dc(mL,g)*dc(mC,g)*dy(j)/(dx(i-1)*dc(mC,g)+dx(i)*dc(mL,g));
    AT( k, g )     = -2*dc(mT,g)*dc(mC,g)*dx(i)/(dy(j+1)*dc(mC,g)+dy(j)*dc(mT,g));
    AC( k, g )     = dx(i)*dy(j)*sr(mC,g) ...
                     - (AL(k-1,g)+AT(k,g)) ...
                     + dc(mC,g)*dy(j)*2*(1-beta1)/(4*(1+beta1)*dc(mC,g)+dx(i)*(1-beta1)) ...
                     + dc(mC,g)*dx(i)*2*(1-beta2)/(4*(1+beta2)*dc(mC,g)+dy(j)*(1-beta2));

    if it == 1
    	Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
    	chi(k,g)   = xi(mC,g);
    else
    	Si( k, g ) = dx(i)*dy(j)*src(cix(i),ciy(j),g); % fixed, group, volume source
    	Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
    	chi(k,g)   = xi(mC,g);
    end
    sct(k,g,:) = dx(i)*dy(j)*sc(mC,g,:);
    siga(k,g)  = dx(i)*dy(j)*ab(mC,g);
    if BCR == 2 % add the right incident source
        Si(k,g) = Si(k,g) + 4*dc(mC,g)*dy(j)*alphaR(j,g)/(4*dc(mC,g)+dx(i));
    end
    if BCB == 2 % add the bottom incident source
        Si(k,g) = Si(k,g) + 4*dc(mC,g)*dx(i)*alphaB(i,g)/(4*dc(mC,g)+dy(j));
    end    

    % TOP RIGHT -----------------------------------------------------------
    i=N; j=M;
    k = i+(j-1)*( N );
    if BCR==1
        beta1=1;
    else
        beta1=0;
    end
    if BCT==1
        beta2=1;
    else
        beta2=0;
    end
    
    % materials
    mL = mt( cix(i-1), ciy(j  ) );
    mB = mt( cix(i  ), ciy(j-1) );
    mC = mt( cix(i  ), ciy(j  ) );

    AL( k-1, g )   = -2*dc(mL,g)*dc(mC,g)*dy(j)/(dx(i-1)*dc(mC,g)+dx(i)*dc(mL,g));
    AB( k-N, g )   = -2*dc(mB,g)*dc(mC,g)*dx(i)/(dy(j-1)*dc(mC,g)+dy(j)*dc(mB,g));
    AC( k, g )     = dx(i)*dy(j)*sr(mC,g) ...
                     - (AL(k-1,g)+AB(k-N,g)) ...
                     + dc(mC,g)*dy(j)*2*(1-beta1)/(4*(1+beta1)*dc(mC,g)+dx(i)*(1-beta1)) ...
                     + dc(mC,g)*dx(i)*2*(1-beta2)/(4*(1+beta2)*dc(mC,g)+dy(j)*(1-beta2));

    if it == 1
    	Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
    	chi(k,g)   = xi(mC,g);
    else
    	Si( k, g ) = dx(i)*dy(j)*src(cix(i),ciy(j),g); % fixed, group, volume source
    	Sh( k, g ) = dx(i)*dy(j)*ns(mC,g);  % put fixed source back
    	chi(k,g)   = xi(mC,g);
    end         
    sct(k,g,:) = dx(i)*dy(j)*sc(mC,g,:);
    siga(k,g)  = dx(i)*dy(j)*ab(mC,g);
    if BCR == 2 % add the right incident source
        Si(k,g) = Si(k,g) + 4*dc(mC,g)*dy(j)*alphaR(j,g)/(4*dc(mC,g)+dx(i));
    end
    if BCT == 2 % add the top incident source
        Si(k,g) = Si(k,g) + 4*dc(mC,g)*dx(i)*alphaT(i,g)/(4*dc(mC,g)+dy(j));
    end    

end % group loop

